﻿/* Copyright 2008 dnAnalytics Project.
 *
 * Contributors to this file:
 * Marcus Cuda
 *
 * Redistribution and use in source and binary forms, with or without modification,
 * are permitted provided that the following conditions are met:
 * 
 * * Redistributions of source code must retain the above copyright notice, this 
 *   list of conditions and the following disclaimer.
 * * Redistributions in binary form must reproduce the above copyright notice, 
 *   this list of conditions and the following disclaimer in the documentation
 *   and/or other materials provided with the distribution.
 * * Neither the name of the dnAnalytics Project nor the names of its contributors
 *   may be used to endorse or promote products derived from this software without
 *   specific prior written permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
 * TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

using System;

namespace dnAnalytics.LinearAlgebra
{
    public partial class HouseholderTransformation
    {
        /// <summary>
        /// Initializes a new instance of the <see cref="HouseholderTransformation"/> class.
        /// </summary>
        /// <param name="vector">The vector to transform.</param>
        public HouseholderTransformation(Vector vector)
        {
            if (vector == null)
            {
                throw new ArgumentNullException("vector");
            }
            mVector = vector.Clone();
            Compute();
        }

        //based on Golub Matrix Computations 3rd edition sections 5.1.3-5.1.4
        private void Compute()
        {
            Vector v = mVector.Clone();
            double sigma = 0;
            for (int i = 1; i < mVector.Count; i++)
            {
                sigma += mVector[i]*mVector[i];
            }

            if (sigma != 0)
            {
                double mu = System.Math.Sqrt(mVector[0]*mVector[0] + sigma);
                if (mVector[0] <= 0)
                {
                    v[0] = mVector[0] - mu;
                }
                else
                {
                    v[0] = -sigma/(mVector[0] + mu);
                }
                double s = (2.0*v[0]*v[0])/(sigma + (v[0]*v[0]));
                for (int i = v.Count - 1; i >= 0; i--)
                {
                    v[i] = v[i]/v[0];
                }
                double w = s*mVector.DotProduct(v);
                mBeta = mVector[0] - v[0]*w;
                mVector = v;
            }
        }
    }
}
